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Abstract 

This  paper  summarizes  the  development  of  a  homogenized  free  energy  model  which  characterizes  the  temper¬ 
ature-dependent  hysteresis  and  constitutive  nonlinearities  inherent  to  relaxor  ferroelectric  materials.  A  kernel 
for  the  model  is  developed  through  mesoscopic  energy  analysis  and  extended  to  provide  macroscopic  constitu¬ 
tive  relations  through  stochastic  homogenization  techniques  based  on  the  assumption  that  certain  underlying 
parameters  are  manifestations  of  underlying  densities  rather  than  constants.  Mechanisms  characterizing  the  de¬ 
crease  in  hysteresis  and  saturation  polarization  polarization  as  temperatures  are  increased  are  constructed  using 
asymptotic  properties  of  the  kernel  which  is  derived  from  statistical  mechanics  tenets.  Attributes  of  the  model 
are  illustrated  through  comparison  with  PMN-PT-BT  data. 
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1.  Introduction 

Transducers  employing  relaxor  ferroelectric  materials  are  increasingly  considered  for  applications  ranging 
from  sonar  transduction  to  high  precision  optics  due  to  their  relatively  large  stroke  capabilities  and  nearly 
anhysteretic  behavior  when  employed  within  their  diffuse  transition  region.  However,  the  price  paid  for  these 
superior  performance  capabilities  is  strong  temperature-dependence  in  the  dielectric  properties  of  the  materials. 
As  illustrated  by  quasistatic  PMN-PT-BT  data  in  Figure  1,  a  decrease  of  40°C  near  the  freezing  temperature 
can  lead  from  anhysteretic  to  hysteretic  E-P  relations  with  a  significant  increase  in  the  saturation  polarization. 
Since  heat  produced  during  transduction  is  proportional  to  the  area  of  the  hysteresis  loop,  materials  operating 
in  hysteretic  regimes  exhibit  increased  temperatures  which  subsequently  reduces  the  hysteresis  —  this  process 
proceeds  until  an  equilibrium  temperature  is  reached.  To  place  the  magnitude  of  these  effects  in  perspective,  we 
note  that  a  PMN-driven  flextensional  sonar  transducer  submersed  in  water  experiences  a  temperature  increase  of 
approximately  40°C  before  equilibrium  is  reached  [9].  Although  less  dramatic,  similar  mechanisms  degrade  the 
accuracy  of  relaxor-based  nano-  or  micropositioners  unless  temperatures  are  carefully  regulated.  This  motivates 
the  development  of  comprehensive  material  models  which  characterize  the  temperature-dependent  dielectric 
behavior  in  a  manner  which  facilitates  subsequent  transducer  design  and  model-based  control  implementation. 

A  number  of  the  initial  models  for  the  relaxor  ferroelectric  compound  lead  magnesium  niobate  (PMN),  and 
derivatives  doped  with  lead  titanate  (PT),  barium  titanate  (BT)  and  strontium  titanate  (ST),  focused  on  the 
quantification  of  the  anhysteretic  behavior  above  the  freezing  temperature  Tf  —  e.g.,  see  [2,  3,  4,  5,  6].  Smith  and 
Horn  subsequently  incorporated  the  hysteretic  low-temperature  behavior  through  the  development  of  domain  wall 
models  having  temperature-dependent  coefficients  [13,  14].  Whereas  these  models  are  very  effective  for  symmetric 
loop  behavior,  they  require  a  priori  knowledge  of  turning  points  to  guarantee  biased  minor  loop  closure  which 
precludes  their  use  in  feedback  control  designs  for  transient  dynamics. 
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Figure  1:  Decreasing  hysteresis  and  saturation  polarization  of  PMN-PT-BT  as  a  function  of  increasing  temper¬ 
ature  (from  [13,  14]). 


In  this  paper,  we  extend  the  homogenized  free  energy  framework,  developed  in  [17]  for  ferroelectric  materials, 
to  incorporate  temperature  effects  inherent  to  relaxor  ferroelectric  compounds.  In  Section  2,  we  summarize 
the  previous  fixed-temperature  model  which  is  comprised  of  an  energy-based  kernel  derived  at  the  mesoscopic 
level  with  macroscopic  effective  parameters  determined  through  stochastic  homogenization  techniques.  Ther¬ 
modynamic  properties  of  the  kernel  are  exploited  in  Section  3  to  construct  mechanisms  that  incorporate  the 
temperature-dependent  behavior  illustrated  in  Figure  1.  Attributes  of  the  model  are  illustrated  in  Section  4 
through  a  comparison  with  quasistatic  PMN-PT-BT  data. 


2.  Isothermal  Hysteresis  Model 

To  provide  the  framework  for  incorporating  temperature-dependent  behavior,  we  summarize  first  the  hys¬ 
teresis  model  developed  in  [17]  for  ferroelectric  materials. 


2.1.  Helmholtz  and  Gibbs  Energies 

As  detailed  in  [17],  combination  of  the  internal  energy  derived  under  the  assumption  that  dipoles  align  either 
with  the  applied  field,  or  diametrically  opposite  to  it,  and  classical  entropy  relations  yields  the  temperature- 
dependent  Helmholtz  energy  relation 
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Here  Eh,Tc  and  Ps  respectively  denote  a  bias  field,  the  Curie  temperature  and  the  saturation  polarization.  It 
can  be  easily  verified  that  (1)  incorporates  double  well  behavior  for  T  <  Tc  and  a  single  potential  well  for 
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V(P)=G(0,P) 


Figure  2:  (a)  Helmholtz  energy  ip  and  Gibbs  energy  G  for  increasing  field  E  (E2  >  E\  >0).  (b)  Dependence  of 
the  local  average  polarization  P  on  the  field  E. 

T  >  Tc  in  accordance  with  the  ferroelectric  to  paraelectric  phase  transition  associated  with  materials  exhibiting 
second-order  phase  transitions. 

Whereas  derived  from  statistical  mechanics  tenets,  the  Helmholtz  relation  (1)  is  often  better  suited  for 
qualitative  rather  than  quantitative  characterization  due  to  the  restricted  assumptions  leading  to  its  formulation. 
To  provide  additional  flexibility  and  efficiency  in  fixed-temperature  regimes,  a  second  choice  for  the  Helmholtz 
functional  is  the  piecewise  quadratic  construct 

r  iv(p+pr)2  ,p<-Pi 

%p{P)  =  |  h(p  -  Pr )2  ,  P>  Pi  (2) 

(  h(Pi  -  Pr)  (g  -  Pr)  ,  \P\  <  Pi 

which  can  be  interpreted  as  being  derived  through  Taylor  expansions  at  each  of  the  equilibria  of  (1).  As  illustrated 
in  Figure  2,  Pj  and  Pr  respectively  denote  the  inflection  point  and  polarization  at  which  the  minimum  of  ip 
occurs. 

Recalling  that  the  electrostatic  energy  is  given  by  Ue  =  p  •  E,  where  p  is  a  dipole  in  an  electric  field  E,  we 
construct  the  Gibbs  energy  relation 

G  =  ip-EP  (3) 

where  ip  is  given  by  (2)  or  (3).  The  behavior  of  G  for  increasing  fields  is  illustrated  in  Figure  2(a). 

2.2.  Local  Polarization  Relation  for  Negligible  Thermal  Activation 

For  operating  regimes  in  which  thermal  activation  is  negligible  and  relaxation  times  are  small  compared 
with  drive  frequencies,  the  local  average  polarization  P  is  determined  directly  the  minimization  of  the  Gibbs 
energy  (3).  For  the  Helmholtz  relation  (1),  enforcement  of  the  necessary  condition  =  0  yields 

P{E)  =  P.tanh  (^f) 

where  a  =  and  a{T )  =  .  The  kernel  resulting  from  the  piecewise  quadratic  relation  (2)  has  the  general 

form 

P(E)  =  -E+PrS  (4) 

77 


3 


where  <5  =  —1  for  negatively  oriented  dipoles  and  <)'  =  1  for  those  with  positive  orientation.  To  specify  S  and 
hence  P  more  specifically  in  terms  of  the  initial  dipole  orientations  and  previous  switches,  we  employ  the  Preisach 
notation  and  take 


r  [P(E-,Ec,Om  ,  r(t)  =  0 

{P(E;  Ec,£)](t)  =  <  f  ~Pr  ,  r(t)  ±  0  and  £(max -r(t))  =  -Ec  (5) 

1  f  +  Pr  ,  r(t)  ±  0  and  E  (max.  r(t))  =  Ec 

where 

r  ,  £(o) <  -sc 

[?(£;  Ec,  0]  (0)  =  i  £  i  /f,  <  £(0)  <  Ec 

If  +Pr  ,  ^(0)  > 

defines  the  initial  states  of  the  kernel.  The  local  coercive  field 

Ec  =  v(Pr  ~  Pi)  (6) 

quantifies  the  field  at  which  the  negative  well  ceases  to  exist  and  hence  a  dipole  switch  occurs.  Finally,  if  tf 
designates  the  final  time  under  consideration,  the  set  of  switching  times  is  given  by 

r(t)  =  {f  £  (0 ,tf]  |  E(t)  =  -Ec  or  E(t)  =  Ec}. 

The  behavior  of  the  piecewise  linear  hysteron  (4)  or  (5)  is  illustrated  in  Figure  2(b),  and  this  is  the  kernel  we 
consider  throughout  the  remainder  of  the  development. 


2.3.  Local  Polarization  Relation  Incorporating  Thermal  Activation 

To  include  the  effects  of  thermal  relaxation,  it  is  necessary  to  balance  the  Gibbs  energy  G  with  the  relative 
thermal  energy  kT /V ,  where  k  and  V  respectively  denote  Boltzmann’s  constant  and  a  reference  volume,  through 
the  Boltzmann  relation 

M(G)  =  Ce~GV'kT. 

The  constant  C  is  specified  to  guarantee  integration  to  unity  when  considered  over  all  admissible  dipole  orien¬ 
tations. 

The  local  average  polarization  in  this  case  is  given  by 

P  =  x+ (P+)  +  x- (P -)  (7) 


where  x+  and  X-  respectively  denote  the  fractions  of  dipoles  having  positive  and  negative  orientations  and 
(P+) ,  (P-)  are  the  associated  average  polarizations.  As  detailed  in  [17],  the  latter  are  quantified  by  the  general 
relations 

/oo  r—Pi 

>  Pf*(G)dP  ,  (P-)  =  J 

Evaluation  of  the  integration  constant  subsequently  yields 

r  Pe~G^vlkTdP  [ 


Pfi(G)  dP. 


-Pi 


pe-G(E,P)V/kTdp 
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>  Pi 


I Pi 


-G(E,P)V/kTdp 


(P-)  = 


-Pi 


-G(E,P)V/kTdp 


The  evolution  of  dipole  fractions  are  quantified  by  the  differential  equations 


x+  =  ~P+- x+  +p-+x_ 

X-  =  - P-+X -  +  p+-X+ 
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which  can  be  simplified  to 


x+  =  -P+-X+  +  p_+(l  -  x+) 


through  the  identity  x+  +  x_  =1.  Here 

[Pl  e-G(E,P)V/kTdp 
1  JPr-e  1 

T(T)  ^  e—G(E,P)V/kT  dp  ’  ^  T(T) 

J  P,  -  f 

respectively  denote  the  likelihoods  that  dipoles  switch  from  positive  to  negative,  and  conversely.  In  (8),  e  is 
a  small  positive  constant  and  T  denotes  the  material-dependent  relaxation  time.  As  detailed  in  [17],  the  local 
average  polarization  relation  (7)  converges  to  the  piecewise  linear  relation  (5)  in  the  limit  kT /V  — >  0  of  negligible 
relative  thermal  energy. 

2.4.  Macroscopic  Polarization  Model 

The  local  polarization  relations  (5)  and  (7)  are  derived  under  the  assumption  of  a  uniform  lattice  and  effective 
field.  Hence  they  provide  adequate  approximations  for  certain  single  crystal  compounds  but  neglect  the  effects  of 
polycrystallinity,  material  and  stress  nonhomogeneities,  and  variable  effective  fields  which  contribute  significantly 
to  the  hysteretic  behavior  of  general  polycrystalline  compounds.  We  incorporate  these  effects  by  considering  the 
local  coercive  and  effective  fields  to  be  manifestations  of  underlying  distributions  rather  than  constants.  This 
provides  low-order  macroscopic  polarization  employing  the  relations  (5)  or  (7)  as  a  kernel. 

It  is  demonstrated  in  [17]  that  distributed  coercive  and  effective  field  behavior  can  be  incorporated  through 
the  general  formulation 


-Pl+e 


-Pi 


-G(E,P)V/kTdp 


p  —  Pi+e 


-G(E  ,P)V /  kT  dp 


(8) 


[P(E)](t)=  /  /  i/i  (Ec)p2{Ee)[P{Ee  +  E-Ec,0\(t)dEedEc  (9) 

J  0  J — oo 

where  P  is  given  by  (5)  or  (7).  The  densities  v\  and  v-2  characterizing  the  coercive  and  effective  field  distributions 
must  satisfy  certain  physical  decay  criteria  but  otherwise  can  be  considered  to  be  general.  Details  regarding  the 
formulation  of  the  model  and  identification  techniques  for  estimating  the  general  densities  can  be  found  in  [10, 12]. 

Alternatively,  one  can  consider  parametric  forms  for  the  densities  which  satisfy  the  constraint  that  the  coercive 
field  is  positive  and  both  densities  exhibit  exponential  decay.  In  this  latter  vein,  we  consider  the  choices 

vAEc)  =  Cie-^/Ec)l2cf 

(10) 

vi{Ee)  =  c2e_B'/b. 

It  is  established  in  [1]  that  if  c  is  small  compared  with  Ec,  the  norm  and  variance  for  the  lognormal  distribution 
have  the  approximate  values 

(Ec)  «  Ec  ,  a  «  2 Ec  c . 

These  relations  can  be  employed  to  obtain  initial  parameter  estimates  which  reflect  properties  of  the  measured 
data.  In  the  next  section,  we  detail  the  modifications  required  to  incorporate  temperature-dependence  in  the 
densities  (10)  and  hence  the  polarization  model  (9). 


3.  Thermal  Hysteresis  Model 

Because  the  Helmholtz  relation  (1)  was  constructed  through  a  balance  of  internal  energy  and  relative  entropy 
over  a  representative  volume,  it  encapsulates  the  diminishing  hysteresis  observed  in  materials  exhibiting  second- 
order  phase  transitions  as  temperatures  are  increased  through  the  Curie  point.  However,  its  applicability  for 
characterizing  general  temperature-dependent  behavior  of  the  type  illustrated  in  Figure  1  should  be  considered 
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as  qualitative  rather  than  quantitative  due  to  the  restricted  assumptions  underlying  its  formulation.  Instead, 
we  use  it  to  predict  extensions  to  the  piecewise  quadratic  functional  (2)  to  incorporate  relevant  temperature- 
dependencies. 


3.2.  Model  Formulation 

Analysis  of  the  data  in  Figure  1  indicates  that  as  temperatures  are  increased  through  the  freezing  temperature 
Tf,  the  local  coercive  field  Ec,  local  remanence  value  Pr  and  inflection  point  Pj  all  decrease.  We  focus  first  on 
the  latter. 

For  the  Helmholtz  energy  relation  (1),  the  inflection  point  Pj  is  determined  by  =  0  where  G  is  given  (2). 
Since 

d2G  _  -Eh  EhT  1 
dP 2  Ps  TCPS  1  -  (P/Ps)2  ’ 

this  yields 

±P/(T)  =  ±Pay/l-T/Tc. 

Because  the  coercive  field  Ec  is  the  point  at  which  dipole  switches  occur,  it  follows  that 


EC(T)  =  ^P/(T)  -  ^arctanh(P/ (T)/Ps) 

=  Eh  y/1  -  T/Tc  -  ^arctanlVl  -  T/T< 

-L  r. 


=  Eh\J\  —  T/Tc  — 


EhT 

Tr 


\/l  -  T/Tc  +  1(1  -  T/Tcf/2 


(11) 


=  Eh{l-T/Tcf/2 


1—  IT 
3  Tr 


Finally,  through  the  relation  (6),  the  local  remanence  polarization 
Pr{T)  =  tMH  +  Pl(T) 


Eh 


(1  -  T/Tc)3/2 


1  - 


1  T 
3  Tc 


+  Ps\A  —  T  /Tc 


(12) 


inherits  temperature-dependence  from  both  Ec  and  Pj. 

We  generalize  the  expansions  in  (11)  and  (12)  to  obtain 


EC(T)  =  Eh(  1  -  T/Tcy 


and 

Pr{T)  =  —(1  -  T/Tcy  +  Ps(l  -  T/Tc)q. 
ri 

The  temperature-dependent  coercive  field  is  employed  in  the  density  v\  of  (10)  and  the  temperature-dependent 
remanence  value  is  incorporated  in  the  linear  kernel  P  given  by  (4)  or  (7). 


3.2.  Model  Implementation 

To  implement  the  model  (9),  it  is  necessary  to  approximate  the  integrals  in  a  manner  which  is  sufficiently 
efficient  to  permit  eventual  real-time  control  implementation.  We  accomplish  this  by  discretizing  the  integrals 
via  a  composite  Gaussian  quadrature  rule  followed  by  development  of  a  highly  efficient  matrix  algorithm  for 
evaluating  the  if-then  constructs  encapsulated  in  the  kernel  definitions  (4)  or  (5).  Details  regarding  these 
algorithms  can  be  found  in  [17]. 
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The  assumption  that  densities  exhibit  exponential  decay  permits  us  to  truncate  the  domain.  A  composite 
Gaussian  quadrature  discretization  on  the  resulting  compact  support  yields  the  approximate  polarization  model 

Ni  Nj 

[P(E,  T)](t)  =  ^^[P(£ej+£;£Ci,a](t)e-EVV  MEH/Ee<X)/2e]' v.w.  (13) 

*= i  j= i 

where  Ee.,ECi  denote  the  abscissas  associated  with  respective  quadrature  formulae  and  Vi,Wj  are  the  respective 
weights  —  e.g.,  see  pages  698-699  of  [18]. 

We  now  provide  details  illustrating  the  evaluation  of  the  kernel  [P{Ee.  +  E,  Ec.  ,&)](£).  The  crux  of  the 
algorithm  entails  the  determination  of  the  active  kernel  branch  indicated  by  (5)  while  avoiding  the  inefficiency 
inherent  to  direct  evaluation  of  if-then  statements.  To  motivate  the  algorithm,  we  suppose  that  we  want  to 
find  the  polarization  corresponding  to  the  effective  field  E  =  E(tk)  +  Ee  at  the  kth  time  step  for  a  fixed  variation 
Ee  and  coercive  field  Ec.  If  the  polarization  at  the  previous  (k  —  l)th  time  step  is  on  the  positive  branch,  then 
P(Ek)  is  on  the  positive  branch  if  E  —  (—Ec)  is  positive  and  on  the  negative  branch  if  E  —  (—Ec)  is  negative. 
If  P(Ek-i)  is  on  the  negative  branch,  then  P(Ek)  is  on  the  negative  branch  if  E  —  Ec  is  negative  and  on  the 
positive  branch  if  E  —  Ec  is  positive.  Hence  the  polarization  at  P(Ek)  depends  on  E  +  Ec  if  P(Ek-i)  is  on  the 
positive  branch  and  E  +  Ec  if  P(Ek-i)  is  on  the  negative  branch.  Thus,  we  are  interested  in  the  quantity 

A  =  sign  (E  +  A  0EC) 


where  Ao  denotes  the  polarization  branch  at  the  previous  time  step  and 

1,  x  >  0 

sign(x)  =  {  0,  x  =  0 

—1,  x  <  0  . 

We  now  use  A  in  an  algorithm  to  compute  P.  We  begin  by  defining  the  matrices 
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-1 


and  weight  vectors 


WT  = 


V1  = 


w\e 


-El  /b 

ei 7  •  •  • 


■  -1 

1  •• 

•  1 

•  -1 

1  •• 

•  1 

-ElN./b] 

wN.  e 

Ni  x  N, 


J  lxNi 


Vie-MEcl/Ec)/2c}2  ...  Vn  e-[ln(ECN,/Ec)/2cf 


1  x  Ni 


At  each  time  step  we  compute  P  =  -Ee  +  PrA  and  use  these  values  to  calculate  the  matrix  A  as  given  by 


A  =  sign  (.Ee  +  A0EC). 


Finally,  we  update 

P  =  —Ee  +  A  Pr 
V 

Aq  =  A 

and  perform  the  matrix-vector  multiplication 

P  =  CVTPW 


to  compute  the  sum. 
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4.  Experimental  Validation 


To  illustrate  the  performance  of  the  temperature-dependent  hysteresis  model,  we  use  it  to  characterize  the 
PMN-PT-BT  behavior  plotted  in  Figure  1.  To  construct  the  model,  it  is  necessary  to  estimate  the  parameters 

C 

q  =  [Eh,Tc,Psr),c,b,  —,p,q}T 
V 

through  a  least  squares  fit  to  the  data.  To  designate  a  least  squares  functional,  we  let  {Ek,  Pk},  k  =  1, . . . ,  Nd, 
denote  the  combined  data  collected  at  four  temperatures,  263,  278,  293,  and  313,  and  define  parameter-dependent 
solutions  to  (13)  by  P{Ek,T\q).  The  least  squares  functional 


Nd 

J(q)  =  J2\P(Ek,T-,  q)-Pk  |2 

k= 1 

was  then  minimized  using  the  Matlab  simplex  algorithm  fmins  .m  to  obtain  the  parameter  values  Tc  =  330.8  K, 
Eh  =  7.5  x  105  Mv/m,  b  =  2.1  x  1011,  Psp  =  6.0275  (C/m2)(MVm/C),  |  =  10"7  C/MVm,  c  =  0.16,  p  =  1.4, 
and  q  =  0.85. 

The  resulting  model  fit  is  compared  with  the  quasistatic  (collected  at  1  Hz)  PMN-PT-BT  data  in  Figure  3.  It 
is  observed  that  the  model  accurately  characterizes  the  diminishing  hysteresis  and  saturation  polarization  values 
measured  in  the  material  as  temperatures  are  increased  within  the  diffuse  transition  region.  Hence  it  provides 
an  initial  framework  quantifying  temperature-dependent  relaxor  ferroelectric  behavior  of  transducer  design  and 
model-based  control  implementation. 


T=263  K 


-1  -0.5  0  0.5  1 


T=278 K 


-1  -0.5  0  0.5  1 


Figure  3:  PMN-PT-BT  data  ( - )  and  homogenized  free  energy  model  ( - ).  Abscissas:  electric  field  (MV/m), 

ordinates:  polarization  (C/m2). 


5.  Concluding  Remarks 


The  model  developed  here  extends  the  framework  in  [17]  to  incorporate  temperature-dependent  hysteresis  and 
nonlinear  constitutive  behavior  inherent  to  relaxor  ferroelectric  compounds.  Specifically,  it  employs  asymptotic 
properties  of  an  underlying  mesoscopic  kernel  constructed  using  mean  field  theory  to  incorporate  mechanisms 
which  provide  the  decrease  in  hysteresis  and  saturation  polarization  measured  in  relaxors  as  temperatures  are 
increased  within  the  diffuse  transition  region. 

The  underlying  framework  provides  a  number  of  properties  which  make  it  advantageous  for  material  char¬ 
acterization,  transducer  design,  and  model-based  control  design.  As  detailed  in  [17],  the  framework  provides 
reversibility  and  guarantees  closure  of  biased  minor  loops.  It  also  provides  an  energy  basis  for  certain  extended 
Preisach  models  [15]  and  establishes  a  unified  framework  for  characterizing  hysteresis  and  constitutive  nonlinear¬ 
ities  in  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  Finally,  the  fact  that  it  can  be  efficiently  inverted 
provides  the  capability  for  constructing  inverse  filters  for  linear  feedback  or  feedforward  control  designs  [7,  8,  11]. 
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